// topology.cpp: implementation of the topology class.
//
//////////////////////////////////////////////////////////////////////

#include <time.h>

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <deque>
#include <map>
#include <algorithm>

#ifdef _UNIX_
#include <sys/time.h>
#else
#include <sys/timeb.h>
#endif

#ifndef CLOCKS_PER_SEC
#define CLOCKS_PER_SEC 1.0e+06
#endif

#define DUPLICATE_EDGE   10

#include "Vec3D.hpp"
#include "topology.hpp"
#include "mgIndexMap.hpp"
#include "topology/blockalloc.h"

#include "topology_struc.hpp"
#include "topologyTree.h"

using namespace FTools;
using namespace std;


//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

// cTopology
//////////////////////////////////////////////////////////////////////
cTopology::cTopology()
{
  NodeVector.clear();
  NumNode = 0;

  ElemVector.clear();
  NumElem = 0;

  EdgeVector.clear();
  NumEdge = 0;
  NumSurfaces = 0;
  feedbackFunction = NULL;
  EdgeList = NULL;

}

// ~cTopology
//////////////////////////////////////////////////////////////////////
cTopology::~cTopology()
{
  Clear ( );
}

// Clean
//////////////////////////////////////////////////////////////////////
void cTopology::Clear ()
{
  NodeVector.clear();
  NumNode = 0;

  ElemVector.clear();
  NumElem = 0;

  EdgeVector.clear();
  NumEdge = 0;

  BoundList.clear ( );

  if (EdgeList != NULL)
    EdgeList->clear ( );

  if (nodeBlock.size() > 0)
    nodeBlock.clear ( );
  if (edgeBlock.size() > 0)
    edgeBlock.clear ( );
  if (elemBlock.size() > 0)
    elemBlock.clear ( );

  NumSurfaces = 0;
  feedbackFunction = NULL;
}


// InsertMesh
//////////////////////////////////////////////////////////////////////
int cTopology::InsertMesh (int n_node, int n_elem, double *coords, int *conn,
                           void (*func) (char *), bool invert_conn /* = false */)
{
  int        i, j, index;
  TopNode_  *currNode;
  TopElem_  *currElem;

  Clear ( );

  feedbackFunction = func;

  NumNode = n_node;
  NumElem = n_elem;

  if (feedbackFunction != NULL)
    feedbackFunction ("Filling node and element vectors...");

  if (n_node == 0 || n_elem == 0) return 0;

  // fill nodes
  for (i = 0; i < n_node; ++i)
  {
    currNode = nodeBlock.get();
    NodeVector.push_back(currNode);
    NodeVector[i]->coords[0] = coords[i*3+0];
    NodeVector[i]->coords[1] = coords[i*3+1];
    NodeVector[i]->coords[2] = coords[i*3+2];
    NodeVector[i]->active    = 1;
  }

  // boudbox
  for (j = 0; j < 3; j++)
    max[j] = min[j] = NodeVector[0]->coords[j];
  for (i = 1; i < n_node; i++)
  {
    for (j = 0; j < 3; j++)
    {
      if (NodeVector[i]->coords[j] < min[j])
        min[j] = NodeVector[i]->coords[j];
      if (NodeVector[i]->coords[j] > max[j])
        max[j] = NodeVector[i]->coords[j];
    }
  }


  // fill elements
  index = 0;
  if (invert_conn)  // invert connectivity 
  {
    for (i = 0; i < n_elem; i++)
    {
      currElem = elemBlock.get ();
      ElemVector.push_back(currElem);
      ElemVector[i]->type = conn[index];
      // for all nodes
      ElemVector[i]->inc[0] = conn[index+1];
      for (j = 1; j < conn[index]; j++)
        ElemVector[i]->inc[j] = conn[index + conn[index] - j + 1];

      ElemVector[i]->active = 1;
      ElemVector[i]->valid  = 1;
      ElemVector[i]->id     = i;
      index = index + conn[index] + 1;
    }
  }
  else
  {
    for (i = 0; i < n_elem; i++)
    {
      currElem = elemBlock.get ();
      ElemVector.push_back(currElem);
      ElemVector[i]->type = conn[index];
      // for all nodes
      for (j = 0; j < conn[index]; j++)
        ElemVector[i]->inc[j] = conn[index+1+j];

      ElemVector[i]->active = 1;
      ElemVector[i]->valid  = 1;
      ElemVector[i]->id     = i;
      index = index + conn[index] + 1;
    }
  }


  if (feedbackFunction != NULL)
    feedbackFunction ("Building and filling edge vector...");
  // Build the list of edges
  BuildEdgeList ( );

  if (feedbackFunction != NULL)
    feedbackFunction ("Building and filling adjacent edges to nodes...");
  // Adjacent Edges to Node
  BuildAdjEdgesToNode ( );

  if (feedbackFunction != NULL)
    feedbackFunction ("Identifing differente surfaces...");
  // Identify different surfaces
  //IdentifySurfaces ( ); 

  BuildOrientation ( );

  if (feedbackFunction != NULL)
    feedbackFunction ("Identifing boundary edges...");
  // Build the boundary edge list
  BuildEdgesBoundary ( );

  // SetValidFaces ( );

  if (feedbackFunction != NULL)
    feedbackFunction ("Computing normals...");
  // Build normals of elements and nodes
  UpdateNormals ( );


  // set local tolerance
  SetLocalToler ( );

  if (feedbackFunction != NULL)
    feedbackFunction (" ");

  return 1;
}


// BoundBox
//////////////////////////////////////////////////////////////////////
void cTopology::BoundBox (double bmin[3], double bmax[3])
{
  int i;
  for (i = 0; i < 3; i++)
  {
    bmin[i] = min[i];
    bmax[i] = max[i];
  }
}


// BuildNormals
//////////////////////////////////////////////////////////////////////
void cTopology::UpdateNormals (void)
{
  int    i, j, *n_adj, id_node;
  Vec3D  u, v, w;

  //normal of elements
  for (i = 0; i < NumElem; ++i)
  {
    if (ElemVector[i]->type < 3)
      continue;

    u = NodeVector[ElemVector[i]->inc[1]]->coords - NodeVector[ElemVector[i]->inc[0]]->coords;
    v = NodeVector[ElemVector[i]->inc[2]]->coords - NodeVector[ElemVector[i]->inc[0]]->coords;
    w = CrossProd (u, v);
    //if (ElemVector[i]->type == 4)
    //{
    //  Vec3D w1, w2;
    //  u = NodeVector[ElemVector[i]->inc[2]]->coords - NodeVector[ElemVector[i]->inc[0]]->coords;
    //  v = NodeVector[ElemVector[i]->inc[3]]->coords - NodeVector[ElemVector[i]->inc[0]]->coords;
    //  w1 = CrossProd (u, v);
    //  w1 = w1.Normalize();
    //  w2 = w.Normalize();
    //  w = (w1 + w2) / 2.0;
    //}
    
    ElemVector[i]->normal = w.Normalize ( );

  }

  /* update normal to points */
  n_adj = new int[NumNode];
  memset (n_adj, 0, NumNode * sizeof (int));
  for (i = 0; i < NumNode; ++i)
  { 
    n_adj[i] = 0;
    NodeVector[i]->normal[0] = 0.0;
    NodeVector[i]->normal[1] = 0.0;
    NodeVector[i]->normal[2] = 0.0;
  }

  for (i = 0; i < NumElem; ++i)
  {
    if (ElemVector[i]->type < 3)
      continue;

    for (j = 0; j < ElemVector[i]->type; ++j)
    {
      id_node = ElemVector[i]->inc[j];
      NodeVector[id_node]->normal += ElemVector[i]->normal;
      n_adj[id_node]++;
    }
  }

  /* average */
  for (i = 0; i < NumNode; i++)
    NodeVector[i]->normal /= n_adj[i];

  delete []n_adj;

  // normal of nodes

}

//////////////////////////////////////////////////////////////////////
int cTopology::GetCoordNode (int id, double coord[3])
{
  if (id > NumNode-1)
    return 0;
  coord[0] = NodeVector[id]->coords[0];
  coord[1] = NodeVector[id]->coords[1];
  coord[2] = NodeVector[id]->coords[2];
  return (1);
}

//////////////////////////////////////////////////////////////////////
void cTopology::SetCoordNode (int id, double *new_coord)
{
  if (id > NumNode-1)
    return;

  NodeVector[id]->coords[0] = new_coord[0];
  NodeVector[id]->coords[1] = new_coord[1];
  NodeVector[id]->coords[2] = new_coord[2];
}

//////////////////////////////////////////////////////////////////////
int  cTopology::IsBoundaryNode (int id)
{
  if (id > NumNode-1)
    return 0;

  return NodeVector[id]->bound;
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetNormalNode (int id, double normal[3])
{
  if (id > NumNode-1)
    return 0;

  normal[0] = NodeVector[id]->normal[0];
  normal[1] = NodeVector[id]->normal[1];
  normal[2] = NodeVector[id]->normal[2];
  return 1;
}

//////////////////////////////////////////////////////////////////////
void cTopology::SetPtsInative (int id)
{
  if (id > NumNode-1)
    return;

  NodeVector[id]->active = 0;
}

//////////////////////////////////////////////////////////////////////
int cTopology::NumAdjNodeNode (int id)
{
  if (id > NumNode-1)
    return 0;
  return ((int) NodeVector[id]->AdjEdges.size());
}

//////////////////////////////////////////////////////////////////////
void cTopology::AdjNodeNode (int id, int *adj_nodes)
{
  if (id > NumNode-1)
    return;

  int size = (int) NodeVector[id]->AdjEdges.size();
  int i;
  for (i = 0; i < size; ++i)
  {
    if (NodeVector[id]->AdjEdges[i]->id[0] != id)
      adj_nodes[i] = NodeVector[id]->AdjEdges[i]->id[0];
    else
      adj_nodes[i] = NodeVector[id]->AdjEdges[i]->id[1];
  }
}

//////////////////////////////////////////////////////////////////////
int cTopology::OrientAdjNodeNode (int id, int *adj_nodes)
{
  if (id > NumNode-1)
    return 0;

  int size = (int) NodeVector[id]->AdjEdges.size();
  int i;

  // current edge
  int curr_edge = -1; 
  if (NodeVector[id]->bound)  // get las 
  {
    for (int j = 0; j < (int) NodeVector[id]->AdjEdges.size(); ++j)
    {
      if (NodeVector[id]->AdjEdges[j]->id[0] == id)
      {
        curr_edge = NodeVector[id]->AdjEdges[j]->i;
        break;
      }
    }
  }
  else  // any edge can be selected
  {
    curr_edge = NodeVector[id]->AdjEdges[0]->i;
  }

  // find adjnodes
  for (i = 0; i < size; ++i)
  {
    if (EdgeVector[curr_edge]->id[0] != id)
      adj_nodes[i] = EdgeVector[curr_edge]->id[0];
    else
      adj_nodes[i] = EdgeVector[curr_edge]->id[1];

    if (i == size-1) // avoid last search
      break;

    curr_edge = getAdjEdgeNextNodeEdge (id, curr_edge);
  }
  

  return 1;
}


//////////////////////////////////////////////////////////////////////
int cTopology::OrientAdjEdgeToNode (int id, int *adj_edges)
{
  if (id > NumNode-1)
  return 0;

  int size = (int) NodeVector[id]->AdjEdges.size();
  int i;

  // current edge
  int curr_edge = -1; 
  if (NodeVector[id]->bound)  // get las 
  {
    TopEdge_ *ptr_edge;
    for (int j = 0; j < (int) NodeVector[id]->AdjEdges.size(); ++j)
    {
      ptr_edge = NodeVector[id]->AdjEdges[j];
      if (ptr_edge->id[0] == id && ptr_edge->type == 1)
      {
        curr_edge = ptr_edge->i;
        break;
      }
    }
  }
  else  // any edge can be selected
  {
    curr_edge = NodeVector[id]->AdjEdges[0]->i;
  }

  // find adjnodes
  for (i = 0; i < size; ++i)
  {
    adj_edges[i] = curr_edge;

    if (i == size-1) // avoid last search
      break;

    curr_edge = getAdjEdgeNextNodeEdge (id, curr_edge);
  }
  

  return 1;
}

//////////////////////////////////////////////////////////////////////
int  cTopology::IsActNode  (int id)
{
  if (id > NumNode-1)
    return 0;

  return (NodeVector[id]->active);
}


// ReturnNodeId 
//////////////////////////////////////////////////////////////////////
int cTopology::ReturnNodeId (double x, double y, double z)
{
  Vec3D pts (x, y, z);

  int    i, idmin = -1;
  double dist, mindist;

  mindist = (NodeVector[0]->coords - pts).Magnitude();
  idmin = 0;

  for (i = 1; i < NumNode; ++i)
  {
    dist = (NodeVector[i]->coords - pts).Magnitude();
    if (dist < mindist)
    {
      mindist = dist;
      idmin = i;
    }
  }
 
 return idmin;
}

//////////////////////////////////////////////////////////////////////
void cTopology::SetNodeInfo (int id, int flag)
{
  if (id > NumNode-1)
    return;
  NodeVector[id]->extraInfo = flag;
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetNodeInfo (int id)
{
  if (id > NumNode-1)
    return 0;
  return (NodeVector[id]->extraInfo);
}


//////////////////////////////////////////////////////////////////////
int cTopology::GetEdge (int i, int *ei, int *ej)
{
  if (i > NumEdge-1)
    return 0;

  *ei = EdgeVector[i]->id[0];
  *ej = EdgeVector[i]->id[1];
  return 1;
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetEdge( int ei, int ej )
{
  if (ei > NumNode-1)
    return -1;
  if (ej > NumNode-1)
    return -1;

  int n_adj_edge = NodeVector[ei]->AdjEdges.size();
  for (int i = 0; i < n_adj_edge; ++i)
  {
    if (NodeVector[ei]->AdjEdges[i]->id[0] == ej || 
        NodeVector[ei]->AdjEdges[i]->id[1] == ej )
        return NodeVector[ei]->AdjEdges[i]->i;
  }
  

  return -1;
}


//////////////////////////////////////////////////////////////////////
bool cTopology::IsBoundaryEdge( int id )
{
  if (id > NumEdge-1)
    return false;
  if (EdgeVector[id]->type == 1)
    return true;
  return false;
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetEdgeSize (int i, double mid[3], double *size)
{
  if (i > NumEdge-1)
    return 0;
  Vec3D  middle (0.0, 0.0, 0.0);
  
  middle = (NodeVector[EdgeVector[i]->id[0]]->coords + NodeVector[EdgeVector[i]->id[1]]->coords) * 0.5;
  mid[0] = middle[0];
  mid[1] = middle[1];
  mid[2] = middle[2];
  *size = EdgeVector[i]->size;
  return 1;
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetBoundEdge (int i, int *ei, int *ej)
{
  if (i > NumEdgeBound-1)
    return -1;

  *ei = BoundVector[i]->id[0];
  *ej = BoundVector[i]->id[1];
  return BoundVector[i]->i;
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetBEdgeSize (int i, double mid[3], double *size)
{
  if (i > NumEdgeBound-1)
    return 0;
  Vec3D  middle (0.0, 0.0, 0.0);
  
  middle = (NodeVector[BoundVector[i]->id[0]]->coords + NodeVector[BoundVector[i]->id[1]]->coords) * 0.5;
  mid[0] = middle[0];
  mid[1] = middle[1];
  mid[2] = middle[2];
  *size = BoundVector[i]->size;
  return 1;
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetIsEdgeValid (int i)
{
  if (i > NumEdge-1)
    return 0;
  return EdgeVector[i]->valid;
}

//////////////////////////////////////////////////////////////////////
double cTopology::GetSmalestEdge (void)
{
  static double small_size = 0.0;
  int j;

  for (j = 0; j < NumEdge; ++j)
  {
    if (EdgeVector[j]->size == 0.0)
      continue;
    if (small_size == 0.0 || EdgeVector[j]->size < small_size)
      small_size = EdgeVector[j]->size;
  }

  return small_size;
}

//////////////////////////////////////////////////////////////////////
void cTopology::SetEdgeInfo (int id, int flag)
{
  if (id > NumEdge-1)
    return;
  EdgeVector[id]->extraInfo = flag;
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetEdgeInfo (int id)
{
  if (id > NumEdge-1)
    return 0;
  return (EdgeVector[id]->extraInfo);
}


//////////////////////////////////////////////////////////////////////
int cTopology::GetElemNNodes (int id)
{
  if (id > NumElem-1)
    return 0;

  return ElemVector[id]->type;
}

//////////////////////////////////////////////////////////////////////
int* cTopology::GetElemConn (int id)
{
  if (id > NumElem-1)
    return 0;

  return ElemVector[id]->inc;
}
 
//////////////////////////////////////////////////////////////////////
int cTopology::GetElemNorm (int id, double normal[3])
{
  if (id > NumElem-1)
    return 0;
  
  normal[0] = ElemVector[id]->normal[0];
  normal[1] = ElemVector[id]->normal[1];
  normal[2] = ElemVector[id]->normal[2];

  return 1;
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetElemNormArea( int id, double normal[3] )
{
  if (id > NumElem-1)
    return 0;
  Vec3D  u, v, w;

  u = NodeVector[ElemVector[id]->inc[1]]->coords - NodeVector[ElemVector[id]->inc[0]]->coords;
  v = NodeVector[ElemVector[id]->inc[2]]->coords - NodeVector[ElemVector[id]->inc[0]]->coords;
  w = CrossProd (u, v);

  normal[0] = w[0];
  normal[1] = w[1];
  normal[2] = w[2];
  return 1;
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetElemIdEdge (int id, int pos_adj)
{
  if (id > NumElem-1)
    return 0;
  if (pos_adj >= ElemVector[id]->type)
    return -1;

  return (ElemVector[id]->AdjEdges[pos_adj]->i);
}
 
//////////////////////////////////////////////////////////////////////
int cTopology::GetIsElemValid (int id)
{
  if (id > NumElem-1)
    return 0;

  return ElemVector[id]->valid;
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetElemCenter (int id, double center[3])
{
  int i;
  if (id > NumElem-1)
    return 0;
  Vec3D  middle (0,0,0);

  for (i = 0; i < ElemVector[id]->type; ++i)
  {
    middle += NodeVector[ElemVector[id]->inc[i]]->coords;
  }
  middle /= ElemVector[id]->type;

  center[0] = middle[0];
  center[1] = middle[1];
  center[2] = middle[2];

  return 1;
}

//////////////////////////////////////////////////////////////////////
int cTopology::IsActElem (int id)
{
  if (id > NumElem-1)
    return 0;

  return ElemVector[id]->active;
}


//////////////////////////////////////////////////////////////////////
void cTopology::SetElemInfo (int id, int flag)
{
  if (id > NumElem-1)
    return;
  ElemVector[id]->extraInfo = flag;
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetElemInfo (int id)
{
  if (id > NumElem-1)
    return 0;
  return (ElemVector[id]->extraInfo);
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetNumLoops (void)
{
  return (int) BoundList.size ();
}

//////////////////////////////////////////////////////////////////////
int cTopology::GetLoop (int id, int *npts, int **idpts)
{
  int i;
  int nloops = (int) BoundList.size ();

  if (id >= nloops)
    return 0;

  int  size = (int) BoundList[id]->Edges.size();

  *npts = size;
  *idpts = new int[size];

  for (i = 0; i < size; ++i)
    (*idpts)[i] = BoundList[id]->Edges[i]->id[0];

  return 1;
}


// NumAdjEdgeToNode
//////////////////////////////////////////////////////////////////////
int cTopology::NumAdjEdgeToNode (int id)
{
  if (id > NumNode-1)
    return 0;

  return (int) NodeVector[id]->AdjEdges.size();
}

// AdjEdgeToNode
//////////////////////////////////////////////////////////////////////
int cTopology::AdjEdgeToNode (int id, int pos)
{
  if (id > NumNode-1)
    return -1;

  if ((int) NodeVector[id]->AdjEdges.size() == 0 ||
      (int) NodeVector[id]->AdjEdges.size() < pos+1 )
    return -1;

  return NodeVector[id]->AdjEdges[pos]->i;
}

// OppAdjEdgeToNode
//////////////////////////////////////////////////////////////////////
int cTopology::OppAdjEdgeToNode (int id, int *ne, int **edges)
{
  int i, j, ei, ej;

  if (id > NumNode-1)
    return 0;

  if ((int) NodeVector[id]->AdjEdges.size() == 0)
    return 0;

  int nelm, *elem;
  vector <int> tmpedges;
  AdjElemToNode (id, &nelm, &elem);
  for (i = 0; i < nelm; ++i)
  {
    int nadjedge = NumAdjEdgeToElem (elem[i]);
    for (j = 0; j < nadjedge; ++j)
    {
      int adjedge = AdjEdgeToElem (elem[i], j);
      AdjNodeToEdge (adjedge, &ei, &ej);
      if (!(ei == id || ej == id))
        tmpedges.push_back (adjedge);
    }
  }
  delete []elem;

  // get adjacent edges
  if (tmpedges.size () == 0)
  {
    *ne = 0;
    *edges = NULL;
    return 0; 
  }

  *ne = (int) tmpedges.size ();
  *edges = new int[*ne];
  for (i = 0; i < *ne; ++i)
    (*edges)[i] = tmpedges[i];

  return 1;
}

// GetAdjNodeToEdge
//////////////////////////////////////////////////////////////////////
int cTopology::AdjNodeToEdge (int id, int *ei, int *ej)
{
  if (id > NumEdge-1)
    return 0;

  *ei = EdgeVector[id]->id[0];
  *ej = EdgeVector[id]->id[1];
  return 1;
}

// AdjElemToNode
//////////////////////////////////////////////////////////////////////
int cTopology::AdjElemToNode (int id, int *ne, int **elem)
{
  if (id > NumNode-1)
    return 0;

  map <int, int> adj_elem;

  // for all adjacent edges
  int i, j, nedges = NumAdjEdgeToNode (id);
  for (i = 0; i < nedges; ++i)
  {
    int adjedge  = AdjEdgeToNode (id, i);
    int nadjelem = NumAdjElemToEdge (adjedge);
    for (j = 0; j < nadjelem; ++j)
    {
      int adjelem = AdjElemToEdge (adjedge, j);
      if (adj_elem.find (adjelem) == adj_elem.end())
        adj_elem.insert (pair <int, int> (adjelem, adjelem));
    }
  }

  // any adj_elems
  if (adj_elem.size() == 0)
  {
    *ne = 0;
    *elem = NULL;
    return 0;
  }

  *ne = (int) adj_elem.size();
  *elem = new int[*ne];
  map <int,int>::iterator ii;
  for (ii = adj_elem.begin(), i = 0; ii != adj_elem.end(); ++i, ++ii)
  {
    (*elem)[i] = (*ii).second;
  }

  return 1;
}


// NumAdjElemToEdge
//////////////////////////////////////////////////////////////////////
int cTopology::NumAdjElemToEdge (int id)
{
  if (id > NumEdge-1)
    return 0;
  return (int) EdgeVector[id]->AdjElms.size();
}

// AdjElemToEdge
//////////////////////////////////////////////////////////////////////
int cTopology::AdjElemToEdge (int id, int pos)
{
  if (id > NumEdge-1)
    return -1;
  int size = (int) EdgeVector[id]->AdjElms.size();
  if (size == 0 || size < pos+1) 
    return -1;

  return EdgeVector[id]->AdjElms[pos]->id;
}

// NumAdjEdgeToElem
//////////////////////////////////////////////////////////////////////
int cTopology::NumAdjEdgeToElem (int id)
{
  if (id > NumElem-1)
    return 0;
  return (int) ElemVector[id]->AdjEdges.size ();
}

// AdjEdgeToElem
//////////////////////////////////////////////////////////////////////
int cTopology::AdjEdgeToElem (int id, int pos)
{
  if (id > NumElem-1)
    return -1;

  if ((int) ElemVector[id]->AdjEdges.size () == 0 ||
      (int) ElemVector[id]->AdjEdges.size () < pos+1 )
   return -1;

  return ElemVector[id]->AdjEdges[pos]->i;
}

// GetVolume 
//////////////////////////////////////////////////////////////////////
double cTopology::GetVolume ( )
{
  int i;
  Vec3D p0, p1, p2;
  double volume = 0.0;

  // loop for all edges in the elements
  for (i = 0; i < NumElem; i++)
  {
    p0 = NodeVector[ElemVector[i]->inc[0]]->coords;
    p1 = NodeVector[ElemVector[i]->inc[1]]->coords;
    p2 = NodeVector[ElemVector[i]->inc[2]]->coords;

    volume = volume + 1.0 / 6.0 * TripleProd (p0, p1, p2);

    if (ElemVector[i]->type == 4)
    {
      p0 = NodeVector[ElemVector[i]->inc[0]]->coords;
      p1 = NodeVector[ElemVector[i]->inc[2]]->coords;
      p2 = NodeVector[ElemVector[i]->inc[3]]->coords;

      volume = volume + 1.0 / 6.0 * TripleProd (p0, p1, p2);
    }

  }

  return volume;
}


// 
//////////////////////////////////////////////////////////////////////
void cTopology::SetNodeInfo (int id, void *info)
{
  if (id > NumNode-1)
    return;
  NodeVector[id]->info = info;
}

// 
//////////////////////////////////////////////////////////////////////
void cTopology::GetNodeInfo (int id, void **info)
{
  if (id > NumNode-1)
    return;
  *info = NodeVector[id]->info;
}

// 
//////////////////////////////////////////////////////////////////////
void cTopology::SetEdgeInfo (int id, void *info)
{
  if (id > NumEdge-1)
   return;
  EdgeVector[id]->info = info;
}

// 
//////////////////////////////////////////////////////////////////////
void cTopology::GetEdgeInfo (int id, void **info)
{
  if (id > NumEdge-1)
    return;
  *info = EdgeVector[id]->info;
}

// 
//////////////////////////////////////////////////////////////////////
void cTopology::SetElemInfo (int id, void *info)
{
  if (id > NumElem-1)
    return;
  ElemVector[id]->info = info;
}

// 
//////////////////////////////////////////////////////////////////////
void cTopology::GetElemInfo (int id, void **info)
{
  if (id > NumElem-1)
    return;
  *info = ElemVector[id]->info;
}



//////////////////////////////////////////////////////////////////////
// Private methods
//////////////////////////////////////////////////////////////////////

#include <set>

//////////////////////////////////////////////////////////////////////////
// local class to store adjacent nodes
class s3dAdjNodes
{
public:
  s3dAdjNodes(void) {};
  ~s3dAdjNodes(void)
  {
    m_adjNodes.clear();
  }

  void insert (int node)
  {
    m_adjNodes.insert(node);
  }

  int size (void)
  {
    return (int) m_adjNodes.size();
  }

  set <int> m_adjNodes; 
};


// BuildEdgeList
//////////////////////////////////////////////////////////////////////
int cTopology::BuildEdgeList (void)
{
  int      id1, id2, id3 = -1;
  int      i, j, curr_step = 0;
  TopEdge_  *currEdge;
  Vec3D    vEdge;
  double   steps = 0.05;
  char     message[126];

  // clear edge list
  if (EdgeList != NULL)
    delete EdgeList;
  EdgeList = new mshsurfIndexMap;
  EdgeList->clear ( );
  EdgeVector.clear();

  // loop for all edges in the elements
  for (i = 0; i < NumElem; i++)
  {
    // print to feedback function
    if (feedbackFunction != NULL)
    {
      if (i*1.0/NumElem > curr_step * steps)
      {
        sprintf (message, "Building and filling edge vector (%3.1f %%)...", curr_step * steps * 100.0);
        feedbackFunction (message);
        ++curr_step;
      }
    }
    
    // for all 
    for (j = 0; j < ElemVector[i]->type; j++)
    {
      id1 = ElemVector[i]->inc[j];
      id2 = ElemVector[i]->inc[(j+1)%ElemVector[i]->type];
    
      currEdge = (TopEdge_ *) EdgeList->find (id1, id2, id3);

      // try to find a edge
      if (currEdge == NULL)  // new edge
      {
        TopEdge_  *newEdge = edgeBlock.get();
        EdgeVector.push_back (newEdge);
        newEdge->i = (int) EdgeVector.size() - 1;

        newEdge->id[0] = id1;
        newEdge->id[1] = id2;
        newEdge->valid = 1;
        newEdge->type  = 0;

        // size
        vEdge = NodeVector[id1]->coords - NodeVector[id2]->coords;
        newEdge->size = vEdge.Magnitude ( );

        EdgeList->insert (id1, id2, id3, (void *)newEdge); 

        currEdge = newEdge;
      }
      // update type 
      currEdge->type = (currEdge->type + 1);
//      cout << "Edge type = " << i << " " << currEdge->type << endl;


      // update adjacent elements
      currEdge->AdjElms.push_back (ElemVector[i]);

      // update adjacent edge
      ElemVector[i]->AdjEdges.push_back (currEdge);
      if (id1 == currEdge->id[0])
        ElemVector[i]->cycles.push_back (0);
      else
        ElemVector[i]->cycles.push_back (1);
    }
  }

  //cpu_time = (clock( ) - cpu_time)/(CLOCKS_PER_SEC*1.0);
  //printf("\t\tCPU time Edge 2 ............. %0.5f (s)\n", (double)cpu_time);

  // fill edge vector
  NumEdge = EdgeList->size ();
  // printf ("Numero de arestas = %d\n", NumEdge);

  EdgeList->clear ( );
  delete EdgeList;
  EdgeList = NULL;

  return 1;
}


// BuildAdjEdgesToNode
//////////////////////////////////////////////////////////////////////
int cTopology::BuildAdjEdgesToNode (void)
{
  int i, id0, id1;

  // clear any adjacent list
  for (i = 0; i < NumNode; ++i)
    NodeVector[i]->AdjEdges.clear();

  for (i = 0; i < NumEdge; ++i)
  {
    id0 = EdgeVector[i]->id[0];
    id1 = EdgeVector[i]->id[1];
    NodeVector[id0]->AdjEdges.push_back (EdgeVector[i]);
    NodeVector[id1]->AdjEdges.push_back (EdgeVector[i]);
    if (EdgeVector[i]->type == 1)
    {
      NodeVector[id0]->bound = 1;
      NodeVector[id1]->bound = 1;
    }
  }

  return 1;
}

// BuildEdgesBoundary
//////////////////////////////////////////////////////////////////////
int cTopology::BuildEdgesBoundary  (void)
{
  int i, j, k;

  BoundList.clear();
  BoundVector.clear();

  // look for a boundary edge
  for (i = 0; i < NumEdge; i++)
  {
    if (EdgeVector[i]->flag == 0  &&  EdgeVector[i]->type == 1)
    {
      // cout << "Uma aresta de contorno" << endl;
      RunAlongEdges (EdgeVector[i]);
    }
  }

  // fill boundary edge vector
  int nloops = (int) BoundList.size ();
  for (i = 0, k = 0; i < nloops; ++i)
  {
    int n = (int) BoundList[i]->Edges.size();
    for (j = 0; j < n; ++j, ++k)
      BoundVector.push_back(BoundList[i]->Edges[j]);
  }

   NumEdgeBound = (int) BoundVector.size();

  return 1;
}


// RunAlongEdges 
//////////////////////////////////////////////////////////////////////
void cTopology::RunAlongEdges (TopEdge_ *initEdge)
{
  int k;
 
  TopEdge_ *currEdge, *adjEdge;
  TopBoundary *newBoundary = new TopBoundary;  // new bounday

  // Add initial edge
  newBoundary->Edges.push_back (initEdge);

  currEdge = initEdge;
  currEdge->flag = 1;

  // obtain the boundary from initial edge
  do
  {
    int id_j, n_edge_adj = 0;

    /* pega o no na posicao de j */
    id_j = currEdge->id[1];

    // count how many adjacent boundary edges to current point
    for (k = 0; k < (int) NodeVector[id_j]->AdjEdges.size(); ++k)
    {
      adjEdge = NodeVector[id_j]->AdjEdges[k];
      if (adjEdge->type == 1 && adjEdge != currEdge && adjEdge->flag == 0)
        ++n_edge_adj;
    }

    // it must find only one edge
    if (n_edge_adj != 1)
      break;

    // update the current edge
    for (k = 0; k < (int) NodeVector[id_j]->AdjEdges.size(); ++k)
    {
      adjEdge = NodeVector[id_j]->AdjEdges[k];
      if (adjEdge->type == 1 && adjEdge != currEdge  && adjEdge->flag == 0)
        currEdge = adjEdge;
    }

    currEdge->flag = 1;

    // add current edge to current boundary
    if (currEdge != initEdge)
      newBoundary->Edges.push_back (currEdge);

  } while (currEdge != initEdge);

  // hold current edge boundary
  BoundList.push_back (newBoundary);
}


// SetLocalToler 
//////////////////////////////////////////////////////////////////////
void cTopology::SetLocalToler (void)
{
  int i;
  LocalTol = 0.0;

  for (i = 0; i < NumEdge; ++i)
  {
    if (EdgeVector[i]->size > 0.0 && (LocalTol == 0.0 || EdgeVector[i]->size < LocalTol))
      LocalTol = EdgeVector[i]->size;
  }

  LocalTol *= 0.01;
}


// IdentifySurfaces 
//////////////////////////////////////////////////////////////////////
void cTopology::IdentifySurfaces ( )
{
  int i;

  NumSurfaces = 0;

  // reset all surfaces in elements
  for (i = 0; i < NumElem; ++i)
    ElemVector[i]->surfId = -1;

  // reset all surfaces in elements
  for (i = 0; i < NumElem; ++i)
  {
    if (ElemVector[i]->surfId == -1)
      SetIdSurface (NumSurfaces, ElemVector[i]);
    ++NumSurfaces;
  }
}

// SetIdSurface 
//////////////////////////////////////////////////////////////////////
void cTopology::SetIdSurface (int idSurf, TopElem_ *initElem)
{
  TopEdge_ *currEdge;
  TopElem_ *currElem;
  deque<TopEdge_ *> dequeEdges;

  int i;

  // reset flags in the edges 
  for (i = 0; i < (int) NumEdge; ++i)
    EdgeVector[i]->flag = 0;

  // add all edges in initial element
  for (i = 0; i < (int) initElem->AdjEdges.size(); ++i)
  {
    if (initElem->AdjEdges[i]->type == 2)
      dequeEdges.push_back (initElem->AdjEdges[i]);
  }

  currElem = initElem;
  currElem->surfId = idSurf;

  // 
  while (!dequeEdges.empty ( ))
  {
    // get the first edge in deque
    currEdge = dequeEdges.front ( );

    // Set flag
    currEdge->flag = 1;

    // set adjacent elements to edge

    // add adjacent edges to deque
    currElem = currEdge->AdjElms[0];
    if (currElem->surfId == -1)
    {
      for (i = 0; i < (int) currElem->AdjEdges.size(); ++i)
      {
        if ((currElem->AdjEdges[i]->type == 2) &&
            (currElem->AdjEdges[i]->flag == 0) )
          dequeEdges.push_back (currElem->AdjEdges[i]);
      }
      currEdge->AdjElms[0]->surfId = idSurf;
    }

    currElem = currEdge->AdjElms[1];
    if (currElem->surfId == -1)
    {
      for (i = 0; i < (int) currElem->AdjEdges.size(); ++i)
      {
        if ((currElem->AdjEdges[i]->type == 2) &&
            (currElem->AdjEdges[i]->flag == 0) )
          dequeEdges.push_back (currElem->AdjEdges[i]);
      }
      currEdge->AdjElms[1]->surfId = idSurf;
    }

    // remove the first edge from deque
    dequeEdges.pop_front ( );
  }
}


// SetValidFaces
//////////////////////////////////////////////////////////////////////
void cTopology::SetValidFaces (void)
{
  int i, j;
  
  // reset all surfaces in elements
  for (i = 0; i < NumElem; ++i)
    ElemVector[i]->valid = 1;

  // set all faces as valid
  for (i = 0; i < NumEdge; ++i)
  {
    if (EdgeVector[i]->AdjElms.size() > 2)
    {
      for (j = 0; j < (int) EdgeVector[i]->AdjElms.size(); ++j)
         EdgeVector[i]->AdjElms[j]->valid = 0;
    }
  }
}


// BuildOrientation 
//////////////////////////////////////////////////////////////////////
void cTopology::BuildOrientation (void)
{
  int      i, k, ok = 1;
  TopElem_ *head_elem, *tail_elem; 

  // set all element as valid
  for (i = 0; i < NumElem; ++i)
    ElemVector[i]->valid = 1;

  // try to any inconsistence in the mesh
  for (i = 0; i < NumEdge; ++i)
  {
    if (EdgeVector[i]->AdjElms.size() > 2)
    {
      EdgeVector[i]->valid = 0;
      for (k = 0; k < (int) EdgeVector[i]->AdjElms.size(); ++k)
        EdgeVector[i]->AdjElms[k]->valid = 0;
    }
    else
      EdgeVector[i]->valid = 1;
  }

  //if (!ok)
  //return;


  /* set all element flags to zero */
  for (i = 0; i < NumElem; ++i)
    ElemVector[i]->flag = 0;

  while (1)
  {
    // look for a element with flag = 0
    for (i = 0; i < NumElem; ++i)
    {
      if (ElemVector[i]->flag == 0)
      {
        head_elem = tail_elem = ElemVector[i];
        break;
      }
    }
    if (i == NumElem)
      break;

    /* from the head element, it orients the all elements */
    do
    {
      head_elem->flag = 1;
    
      for (i = 0; i < head_elem->type; i++)
        OrientAdjElement (head_elem->AdjEdges[i], head_elem->cycles[i], &tail_elem);
    
      head_elem = head_elem->next;
    
    } while (head_elem != NULL);

  }


}


/* OrientAdjElement
////////////////////////////////////////////////////////////////////*/
void cTopology::OrientAdjElement (TopEdge_ *edge, int ciclo, TopElem_ **tail_elem)
{
  TopElem_  *curr_elem;
  int       i, pos = -1, n;
  TopEdge_  *tmp_edge;

  if (edge->type == 1) /* boundary */
    return;

  /* both adjacent elements already were passed */
  if (edge->AdjElms[0]->flag && edge->AdjElms[1]->flag)
    return;

  /* get the element to be checked */
  if (edge->AdjElms[0]->flag == 0)
    curr_elem = edge->AdjElms[0];
  else
    curr_elem = edge->AdjElms[1];
  n = curr_elem->type;

  /* find edge position */
  for (i = 0; i < n; i++)
  {
    if (edge == curr_elem->AdjEdges[i])
    {
      pos = i;
      break;
    }
  }

  /* check ciclo */
  if (curr_elem->cycles[pos] == ciclo)
  {
    int tmp;
    /* invert orientation */
    for (i = 0; i < (n-1)/2; i++)
    {
      /* incidency */
      tmp                   = curr_elem->inc[i+1];
      curr_elem->inc[i+1]   = curr_elem->inc[n-i-1];
      curr_elem->inc[n-i-1] = tmp;
    }

     /* invert edge orientation */
    for (i = 0; i < n/2; i++)
    {
#if 1
      /* ciclo */ 
      tmp                      = curr_elem->cycles[i];
      curr_elem->cycles[i]     = curr_elem->cycles[n-i-1];
      curr_elem->cycles[n-i-1] = tmp;
      /* edge */ 
      tmp_edge                   = curr_elem->AdjEdges[i];
      curr_elem->AdjEdges[i]     = curr_elem->AdjEdges[n-i-1];
      curr_elem->AdjEdges[n-i-1] = tmp_edge;
#endif
    }

    for (i = 0; i < n; i++)
    {
      curr_elem->cycles[i] = !curr_elem->cycles[i]; /* invert ciclo */
      if (curr_elem->AdjEdges[i]->type == 1) /* invert id boundary */
      {
        tmp = curr_elem->AdjEdges[i]->id[0];
        curr_elem->AdjEdges[i]->id[0] = curr_elem->AdjEdges[i]->id[1];
        curr_elem->AdjEdges[i]->id[1] = tmp;
      }
    }
  }

  curr_elem->flag = 1;  /* set flag */
  /* orient adj elements */
  (*tail_elem)->next = curr_elem;
  (*tail_elem) = curr_elem;

}

//////////////////////////////////////////////////////////////////////////
int cTopology::getAdjEdgeNextNodeEdge (int node, int edge)
{
  int n_adj_elem = (int) EdgeVector[edge]->AdjElms.size();

  if (n_adj_elem > 2)
    return -1;

  // next node adjacent to node
  int next_node;
  if (EdgeVector[edge]->id[0] == node)
    next_node = EdgeVector[edge]->id[1];
  else
    next_node = EdgeVector[edge]->id[0];


  // find element

  for (int i = 0; i < n_adj_elem; ++i)
  {
    int adj_elem = EdgeVector[edge]->AdjElms[i]->id;
    // for all edges
    int n_edges = (int) ElemVector[adj_elem]->AdjEdges.size();
    for (int j = 0; j < n_edges; ++j)
    {
      int prev = (j+n_edges-1)%n_edges;
      int next = (j+1)%n_edges;
      int curr = ElemVector[adj_elem]->inc[j];
      int curr_next = ElemVector[adj_elem]->inc[next];
      int curr_prev = ElemVector[adj_elem]->inc[prev];


      // find by the corner node and next node
      if (curr == node &&  curr_next == next_node)
        return (ElemVector[adj_elem]->AdjEdges[prev]->i);

      curr_prev = -1;
    }
  }
  return -1;
}


//////////////////////////////////////////////////////////////////////////
int cTopology::writeNF( char *filename )
{
  int index =0, j, i;
  FILE *arquivo;
  int numq4 = 0, numt3 = 0, numq8 = 0, numt6 = 0;

  arquivo = fopen(filename,"w");
  if (arquivo == NULL)
    return 0;

  fprintf(arquivo,"\n%%HEADER\nFile created by program 'Quebra2D' at - AMiranda\n");
  fprintf(arquivo,"\n%%HEADER.VERSION\n1.00\n");
  fprintf(arquivo,"\n%%HEADER.TITLE\n' untitled'\n");

  fprintf(arquivo,"%%NODE\n");
  fprintf(arquivo,"%d\n\n",NumNode);
  fprintf(arquivo,"%%NODE.COORD\n");
  fprintf(arquivo,"%d\n\n",NumNode);

  for(i=0; i<NumNode; i++)
  {
    fprintf(arquivo,"%d    %f    %f    %f\n",i+1, NodeVector[i]->coords[0], 
            NodeVector[i]->coords[1], NodeVector[i]->coords[2]);
  }

  fprintf(arquivo,"\n%%MATERIAL\n1\n");
  fprintf(arquivo,"\n%%MATERIAL.LABEL\n1\n1\t'mat1'\n");          
  fprintf(arquivo,"\n%%MATERIAL.ISOTROPIC\n1\n1\t10000\t0.2\n");  
  fprintf(arquivo,"\n%%THICKNESS\n1\n1\t1\n");                    
  fprintf(arquivo,"\n%%INTEGRATION.ORDER\n1\n1\t2\t2\t2\t2\t2\t2\n");


  index = 0;
  for(i=0; i<NumElem; i++)
  {
    if (ElemVector[i]->type == 3)
      numt3++;
    if (ElemVector[i]->type == 4)
      numq4++;
    if (ElemVector[i]->type == 6)
      numt6++;
    if (ElemVector[i]->type == 8)
      numq8++;

  }


  fprintf(arquivo,"%%ELEMENT\n");
  fprintf(arquivo,"%d\n\n",NumElem);

  index = 0;
  if (numt3 > 0)
  {
    fprintf(arquivo,"%%ELEMENT.T3\n");
    fprintf(arquivo,"%d\n", NumElem);
    for(i=0; i<NumElem; i++)
    {
      if (ElemVector[i]->type == 3)
      {
        fprintf(arquivo,"%d  1  1  1",i+1);
        for(j=0; j<ElemVector[i]->type; j++)
          fprintf(arquivo,"  %d",ElemVector[i]->inc[j] +1);
      }
      fprintf(arquivo,"\n");
    }
  }

  index = 0;
  if (numt6 > 0)
  {
    fprintf(arquivo,"%%ELEMENT.T6\n");
    fprintf(arquivo,"%d\n", numt6);
    for(i=0; i<NumElem; i++)
    {
      if (ElemVector[i]->type == 6)
      {
        fprintf(arquivo,"%d  1  1  1",i+1);
        for(j=0; j<ElemVector[i]->type; j++)
          fprintf(arquivo,"  %d",ElemVector[i]->inc[j]+1);
      }
      fprintf(arquivo,"\n");
    }
  }

  index = 0;
  if (numq4 > 0)
  {
    fprintf(arquivo,"%%ELEMENT.Q4\n");
    fprintf(arquivo,"%d\n", numq4);
    for(i=0; i<NumElem; i++)
    {
      if (ElemVector[i]->type == 4)
      {
        fprintf(arquivo,"%d  1  1  1",i+1);
        for(j=0; j<ElemVector[i]->type; j++)
          fprintf(arquivo,"  %d",ElemVector[i]->inc[j]+1);
      }
      fprintf(arquivo,"\n");
    }
  }

  index = 0;
  if (numq8 > 0)
  {
    fprintf(arquivo,"%%ELEMENT.Q8\n");
    fprintf(arquivo,"%d\n", numq8);
    for(i=0; i<NumElem; i++)
    {
      if (ElemVector[i]->type == 8)
      {
        fprintf(arquivo,"%d  1  1  1",i+1);
        for(j=0; j<ElemVector[i]->type; j++)
          fprintf(arquivo,"  %d",ElemVector[i]->inc[j]+1);
      }
      fprintf(arquivo,"\n");
    }
  }

  fprintf(arquivo,"\n%%END\n");
  fclose(arquivo);

  return 1;
}

//
//////////////////////////////////////////////////////////////////////////
void cTopology::GetNodeVector( double **node_vec )
{
  *node_vec = new double[NumNode*3];
  for (int i = 0; i < NumNode; ++i)
  {
    (*node_vec)[i*3+0] = NodeVector[i]->coords[0];
    (*node_vec)[i*3+1] = NodeVector[i]->coords[1];
    (*node_vec)[i*3+2] = NodeVector[i]->coords[2];
  }
}

//
//////////////////////////////////////////////////////////////////////////
int cTopology::getIdEdgeToNodes( int i_n, int j_n )
{
  // check validity
  if (i_n >= NumNode || j_n >= NumNode)
    return -1;

  // try to find the edge
  int i, j, ni, nj;
  ni = (int) NodeVector[i_n]->AdjEdges.size();
  nj = (int) NodeVector[j_n]->AdjEdges.size();
  for (i = 0; i < ni; ++i)
  {
    for (j = 0; j < nj; ++j)
    {
      if (NodeVector[i_n]->AdjEdges[i]->i == NodeVector[j_n]->AdjEdges[j]->i)
        return NodeVector[i_n]->AdjEdges[i]->i;
    }
  }

  return -1;
}

//
//////////////////////////////////////////////////////////////////////////
bool cTopology::duplicateEdges(int nedges, int *id_edges, bool inter_elem /*= false*/ )
{
  int i;
  vector <int> possible_nodes;
  map    <int, bool> order_nodes;  // id node and flag => true if node not have to be duplicated  
  map    <int, bool>::iterator itr;  
  pair   <map <int, bool>::iterator, bool> status;

  // check ids of edges and if they has two adjacent elements
  for (i = 0; i < nedges; ++i)
  {
    if (id_edges[i] >= NumEdge)
      return false;
    if (EdgeVector[id_edges[i]]->AdjElms.size() != 2)
      return false;

    // create a structure of all node that can be duplicated
    status = order_nodes.insert (make_pair(EdgeVector[id_edges[i]]->id[0], false));
    itr = status.first;
    itr->second = status.second;
    status = order_nodes.insert (make_pair(EdgeVector[id_edges[i]]->id[1], false));
    itr = status.first;
    itr->second = status.second;
  }

  // sort possible nodes
  sort (possible_nodes.begin(), possible_nodes.end());
  
  // select duplicated nodes

  return true;
}

//
//////////////////////////////////////////////////////////////////////////
bool cTopology::addElements( int n_node, int n_elem, double *coords, int *conn, 
                             void (*func) (char *), topologyTree *node_edge_tree,
                             bool invert_conn)
{
  int i, j;
  TopNode_  *currNode;
  TopElem_  *currElem;


  if (node_edge_tree == NULL)
    return false;

  if (n_node == 0 || n_elem == 0)
    return false;

  // create a temporary topology
  cTopology *tmp_top = new cTopology;
  tmp_top->InsertMesh (n_node, n_elem, coords, conn, func, invert_conn);
  // tmp_top->writeNF ("c:\\temp\\tmp_top.dat");

  // Set all nodes and edges flag to -1. This flags will be used to pointer
  // to correspondent node and edge to this->topology.
  // For all boundary edges, find the correspondent id edges in this current
  // top and and set the flag. Also set node flags.
  int new_nodes = 0;
  markEdgesNodesWhenAddElems (tmp_top, node_edge_tree, &new_nodes);

  // add new nodes
  if (new_nodes > 0)
  {
    int id; double coord_[3];

    // update new nodes
    for (i = 0; i < n_node; ++i)
    {
      id = tmp_top->GetNodeInfo(i);
      if (id >= NumNode)
      {
        // update coordinate
        tmp_top->GetCoordNode(i, coord_);
        currNode = nodeBlock.get();
        NodeVector.push_back(currNode);
        NodeVector[id]->coords[0] = coord_[0];
        NodeVector[id]->coords[1] = coord_[1];
        NodeVector[id]->coords[2] = coord_[2];
        NodeVector[id]->active    = 1;

        // update bound box
        for (j = 0; j < 3; j++)
        {
          if (NodeVector[id]->coords[j] < min[j])
            min[j] = NodeVector[id]->coords[j];
          if (NodeVector[id]->coords[j] > max[j])
            max[j] = NodeVector[id]->coords[j];
        }
      }
    }

    NumNode += new_nodes;
  }


  // add new elements
  int n_node_elem, *conn_;
  for (i = 0; i < n_elem; i++)
  {
    n_node_elem = tmp_top->GetElemNNodes (i); 
    conn_ = tmp_top->GetElemConn (i);

    currElem = elemBlock.get ();
    ElemVector.push_back(currElem);
    ElemVector[NumElem+i]->type = n_node_elem;
    for (j = 0; j < n_node_elem; ++j)
      ElemVector[NumElem+i]->inc[j] = tmp_top->GetNodeInfo(conn_[j]);

    ElemVector[NumElem+i]->active = 1;
    ElemVector[NumElem+i]->valid  = 1;
    ElemVector[NumElem+i]->id     = NumElem+i;
  }
  NumElem += n_elem;

  
  // update Edges
  BuildEdgeList ();
  
  BuildAdjEdgesToNode ( );

  BuildOrientation ( );

  BuildEdgesBoundary ( );


  return true;
}

//
//////////////////////////////////////////////////////////////////////////
void cTopology::
markEdgesNodesWhenAddElems( cTopology *tmp_top, topologyTree *node_edge_tree, int *new_nodes)
{
  int i, j, k, m, id, ei, ej, e_i, e_j;
  double coord[3], size;
  vector <int> found, found_ni, found_nj;
  bool   f_i, f_j;

  if (node_edge_tree == NULL || tmp_top == NULL)
    return;

  int n_node = tmp_top->NumNodes();

  // Set all nodes and edges flag to -1. This flags will be used to pointer
  // to correspondent node and edge to this->topology.
  // For all boundary edges, find the correspondent id edges in this current
  // top and and set the flag. Also set node flags.
  for (i = 0; i < n_node; ++i)
    tmp_top->SetNodeInfo(i, -1);
  int nbound = tmp_top->NumBoundEdges ();
  for (i = 0; i < nbound; ++i)
  {
    id = tmp_top->GetBoundEdge(i, &ei, &ej);
    tmp_top->SetEdgeInfo (id, -1);

    // find edges
    tmp_top->GetEdgeSize (id, coord, &size);
    found = node_edge_tree->findEdges (coord[0], coord[1], coord[2],  
      coord[0], coord[1], coord[2]);

    // find nodes
    tmp_top->GetCoordNode (ei, coord);
    found_ni = node_edge_tree->findNodes (coord[0], coord[1], coord[2],  
      coord[0], coord[1], coord[2]);
    tmp_top->GetCoordNode (ej, coord);
    found_nj = node_edge_tree->findNodes (coord[0], coord[1], coord[2],  
      coord[0], coord[1], coord[2]);

    // printf ("***********************************\n");

    // for all found edges
    for (j = 0; j < (int) found.size(); ++j)
    {
      // search for nodes
      e_i = EdgeVector[found[j]]->id[0];
      e_j = EdgeVector[found[j]]->id[1];
      f_i = f_j = false;


      // node i
      for (k = 0; k < (int) found_ni.size(); ++k)
      {
        if (found_ni[k] == e_i || found_ni[k] == e_j)
        {
          // printf ("Node i = %d\n", found_ni[k]+1);
          f_i = true;
          break;
        }
      }

      // node j
      for (m = 0; m < (int) found_nj.size(); ++m)
      {
        if (found_nj[m] == e_i || found_nj[m] == e_j)
        {
          // printf ("Node j = %d\n", found_nj[m]+1);
          f_j = true;
          break;
        }
      }

      // check for searching to set flag to node
      if (f_i == true && f_j == true)
      {
        tmp_top->SetEdgeInfo (id, EdgeVector[found[j]]->i);
       if (found_ni[k] == e_i)
        {
          tmp_top->SetNodeInfo(ei, e_i);
          tmp_top->SetNodeInfo(ej, e_j);
          // printf ("Edge = %d - %d\n", e_i+1, e_j+1);
        }
        else
        {
          tmp_top->SetNodeInfo(ei, e_j);
          tmp_top->SetNodeInfo(ej, e_i);
          // printf ("Edge = %d - %d\n", e_j+1, e_i+1);
        }
        break;
      }  // if
    } // for all found edges
  } // for bound edges

  // count number of nodes and add new ids
  *new_nodes = 0;
  for (i = 0; i < n_node; ++i)
  {
    if (tmp_top->GetNodeInfo(i) == -1)
    {
      tmp_top->SetNodeInfo(i, NumNode + (*new_nodes));
      (*new_nodes)++;
    }
  }

  // printf ("Add new nodes = %d\n", *new_nodes);

}

//
//////////////////////////////////////////////////////////////////////////
bool cTopology::duplicateNodes( int nnodes, int *nodes, bool inter_elem /*= false*/ )
{
  int i, j, n;
  bool duplicateEnd[2] = {false, false};

  // check if end nodes must be duplicated
  int endNodes[2] = {nodes[0], nodes[nnodes-1]};
  for (i = 0; i < 2; ++i)
  {
    // for all adjacent edges to node
    n = (int) NodeVector[endNodes[i]]->AdjEdges.size();
    for (j = 0; j < n; ++j)
    {
      if (NodeVector[endNodes[i]]->AdjEdges[j]->type == 1) // boundary node
      {
        duplicateEnd[i] = true;
        break;
      }
    }
  }

  // get edges from nodes and set as future duplication
  int *edgesToDuplicate = new int[nnodes-1];
  int *oldEdgeInfo      = new int[nnodes-1];
  for (i = 0; i < (nnodes-1); ++i)
  {
    edgesToDuplicate[i] = this->getIdEdgeToNodes(nodes[i], nodes[i+1]);
    oldEdgeInfo[i]      = this->GetEdgeInfo (edgesToDuplicate[i]);
    this->SetEdgeInfo(edgesToDuplicate[i], DUPLICATE_EDGE);
  }



  // release
  for (i = 0; i < (nnodes-1); ++i)
    this->SetEdgeInfo(edgesToDuplicate[i], oldEdgeInfo[i]);
  delete []edgesToDuplicate;
  delete []oldEdgeInfo;

  return true;
}

//
//////////////////////////////////////////////////////////////////////////
void cTopology::duplicateNode( int node )
{
  
  // get number of edges adjacent do node
  int n_adj = (int) NodeVector[node]->AdjEdges.size();

  if (n_adj == 2)
  {

    return;
  }
}

//
//////////////////////////////////////////////////////////////////////////
void cTopology::duplicateNode2Adj( int node )
{


}

//
//////////////////////////////////////////////////////////////////////////
int cTopology::writeMatlab( char *filename, char *open_type)
{
  int index =0, i;
  FILE *f;
  int num_gen = 0;

  f = fopen(filename, open_type);
  if (f == NULL)
    return 0;

  fprintf (f, "X = [ \n");
  for (i = 0; i < NumNode; ++i)
    fprintf (f, "     %10.5f\n", NodeVector[i]->coords[0]);
  fprintf (f, "] \n\n");

  fprintf (f, "Y = [ \n");
  for (i = 0; i < NumNode; ++i)
    fprintf (f, "     %10.5f\n", NodeVector[i]->coords[1]);
  fprintf (f, "] \n\n");

  fprintf (f, "Z = [ \n");
  for (i = 0; i < NumNode; ++i)
    fprintf (f, "     %10.5f\n", NodeVector[i]->coords[2]);
  fprintf (f, "] \n\n");

  fprintf (f, "C = [ \n");


  index = 0;
  for(i=0; i<NumElem; i++)
  {
    if (ElemVector[i]->type == 3)
    {
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[0]+1, ElemVector[i]->inc[1]+1, ElemVector[i]->inc[2]+1);
    }
    else if (ElemVector[i]->type == 4)
    {
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[0]+1, ElemVector[i]->inc[1]+1, ElemVector[i]->inc[2]+1);
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[0]+1, ElemVector[i]->inc[2]+1, ElemVector[i]->inc[3]+1);
    }
    else if (ElemVector[i]->type == 6)
    {
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[0]+1, ElemVector[i]->inc[1]+1, ElemVector[i]->inc[5]+1);
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[1]+1, ElemVector[i]->inc[2]+1, ElemVector[i]->inc[3]+1);
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[3]+1, ElemVector[i]->inc[4]+1, ElemVector[i]->inc[5]+1);
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[1]+1, ElemVector[i]->inc[3]+1, ElemVector[i]->inc[5]+1);
    }
    else if (ElemVector[i]->type == 8)
    {
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[0]+1, ElemVector[i]->inc[1]+1, ElemVector[i]->inc[7]+1);
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[1]+1, ElemVector[i]->inc[2]+1, ElemVector[i]->inc[3]+1);
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[3]+1, ElemVector[i]->inc[4]+1, ElemVector[i]->inc[5]+1);
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[5]+1, ElemVector[i]->inc[6]+1, ElemVector[i]->inc[7]+1);
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[1]+1, ElemVector[i]->inc[3]+1, ElemVector[i]->inc[7]+1);
      fprintf (f, "     %5d %5d %5d; \n", ElemVector[i]->inc[3]+1, ElemVector[i]->inc[5]+1, ElemVector[i]->inc[7]+1);
    }
    else
      num_gen++;
  }

  fprintf (f, "] \n\n");

  fprintf (f, "trisurf (C, X, Y, Z)\n");

  fprintf (f, "hold on \n\n");

  fclose(f);

  return 1;
}

//
//////////////////////////////////////////////////////////////////////////
int cTopology::WriteVTK( char *filename )
{
  int j, i;
  FILE *arquivo = NULL;

  arquivo = fopen(filename,"w");
  if (arquivo == NULL)
    return 0;

  fprintf(arquivo,"# vtk DataFile Version 2.0\n");
  fprintf(arquivo,"Modelos com template 3D\n");
  fprintf(arquivo,"ASCII\n\n");

  fprintf(arquivo,"DATASET UNSTRUCTURED_GRID\n");
  fprintf(arquivo,"POINTS %d float\n",NumNode);

  for(i=0; i<NumNode; i++)
  {
    fprintf(arquivo,"%f    %f    %f\n", NodeVector[i]->coords[0], 
      NodeVector[i]->coords[1], NodeVector[i]->coords[2]);
  }



  fprintf(arquivo,"CELLS %d %d\n", NumElem, NumElem*5);
  for(i=0; i<NumElem; i++)
  {
    fprintf(arquivo,"%d", ElemVector[i]->type);
    for(j=0; j< ElemVector[i]->type; j++)
      fprintf(arquivo,"  %d", ElemVector[i]->inc[j]);
    fprintf(arquivo,"\n");
  }

  fprintf(arquivo,"CELL_TYPES %d\n", NumElem);
  for(i=0; i<NumElem; i++)
  {
    if (ElemVector[i]->type == 3)
      fprintf(arquivo,"5\n");
    else if (ElemVector[i]->type == 4)
      fprintf(arquivo,"9\n");
  }

  fclose(arquivo);

  return 1;
}